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Abstract 

The multigrid one shot method for optimal control problems, governed by elliptic 
systems, is introduced for the infinite dimensional control space. In this case the con- 
trol variable is a function whose discrete representation involves increasing number of 
variables with grid refinement. The minimization algorithm uses Lagrange multipliers 
to calculate sensitivity gradients. A preconditioned gradient descent algorithm is ac- 
celerated by a set of coarse grids. It optimizes for different scales in the representation 
of the control variable on different discretization levels. An analysis which reduces 
the problem to the boundary is introduced. It is used to approximate the two level 
asymptotic convergence rate, to determine the amplitude of the minimization step, 
and the choice of a high pass filter to be used when necessary. The effectiveness of the 
method is demonstrated on a series of test problems. The new method enables the 
solutions of optimal control problems at the same cost of solving the corresponding 
analysis problems just a few times. 
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1 Introduction 


Numerous computational methods have been developed for predicting the performance of 
physical systems. For engineering design purposes a modification of the system configura- 
tion that results in optimal performance is required. However, computations of large scale 
optimal control problems are extremely time consuming and in many cases not practical. 
The effort to overcome such computational difficulties is done in the direction of devel- 
oping faster computers, on the one hand, and in improving the performance of existing 
algorithms, on the other hand. 

An important and difficult class of optimal control problems are optimal shape design 
(OSD) of aerodynamic systems [2, 8, 9, 10, 11]; for example, the design of a wing shape. 
Under certain assumptions, OSD problems can be reduced to simpler optimal boundary 
control (OBC) problems using the small disturbance approximation [6]. The resulting 
problems involve a fixed physical domain with control variables that are defined as boundary 
data. The problem is to minimize a cost function under certain constrains, which are a set 
of PDEs called “state equations”. The cost function is defined to measure the performance 
of the physical system. 

A standard solution process involves an iterative algorithm where each iteration is com- 
posed of two steps. First, the control variables are updated with the “sensitivity gradients” 
which are the gradients of the cost function with respect to the control variables. Then the 
state variables are updated by solving the constraint equation with the new values of the 
control variables. The repeated solution of the constraint PDE makes this computation 
extremely time consuming and in some cases not practical. 

Several methods were developed to calculate the sensitivity gradients. Among them 
is the “adjoint method” [5, 9, 11]. In the adjoint formulation, a Lagrangian is defined 
together with Lagrange multipliers, which are also called “costate variables”. Costate and 
design equations are derived from the variation of the Lagrangian, and together with the 
state equation form the necessary conditions for a minimum. The sensitivity gradients are 
the design equation residuals calculated with the solutions of the state and costate PDEs. 
The adjoint method was first applied to aerodynamic design by A. Jameson in 1988 [9]. A 
multigrid (MG) solver was used to accelerate the convergence of the solution of the state 
and costate equations. This reduces the computational cost of each optimization step to 
O(N) operations (where N is the number of computational grid points) but does not reduce 
the number of iterations to reach the minimum. 

Originally MG methods were developed to accelerate the convergence rate of the nu- 
merical solutions of PDEs [3, 7, 14]. A. Brandt suggested in 1984 [4, page 119] to apply MG 
methods for optimization problems in the framework of a full multigrid (FMG) algorithm 
where the optimization problem should be solved on coarse levels and interpolated to finer 
levels until the finest level is reached. It is further suggested in [4] to treat the optimization 
problem on all levels where on the finer grids the optimization step should be done locally 
if possible and the smooth corrections of the error should be done during the coarse grid 
correction. The whole problem should be solved in one application of the FMG solver. 
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The adjoint and MG methods were combined to solve an optimal control problem, in 
“one shot”, for the finite dimensional control by S. Ta’asan in 1991 [12]. In the finite dimen- 
sional approach, the control variables are represented as a finite sum of some preassigned 
base functions. The main idea in [12] is to represent the state, the costate and the design 
equations on coarser grids with the full approximation scheme ( FA S ) [3] . It is shown in [12] 
that in general the use of Lagrange multipliers is essential to achieve acceleration of the fine 
grid solution process by coarser grids. The algorithm optimizes the control variables on 
coarse grids, and thus, eliminates the repeated solution of fine grid equations in every op- 
timization step. The one shot algorithm was applied successfully to the small disturbance 
approximation of an aerodynamic wing design problem in a subsonic flow by S. Ta.'asan, 
G. Kuruvila and M. D. Salas (1992) [13]. The performance of the algorithm in [13] was a 
substantial improvement in terms of computational cost. However, the performance of the 
finite dimensional one shot algorithm depends on on the choice of base functions and on 
the level on which the different control variables were optimized. 

In this paper we extend the multigrid one shot algorithm to the infinite dimensional 
control. We introduce an analysis which reduces the problem to the boundary. The analysis 
is used to to determine a minimization step which reduces mainly the high frequency errors 
in the control variables. In elliptic systems such a minimization step requires an update of 
the state and costate solutions only in a local area neighboring the boundary. Based on 
the above, two level analysis is done to approximate the convergence performance of the 
algorithm for a given problem and discretization. Computational demonstrations of the 
algorithm are given for a set of test problems in which the PDE constraint is elliptic. 


2 A Single Grid Algorithm for the Solution of Opti- 
mal control Problems 

2.1 Problem Definition 

Let Yl be a bounded open set of IR d with smooth boundary Y and let (ft be a real valued 
function on fi. Let U and W be Hilbert spaces of real valued functions which are defined 
on T and Yl respectively. 

The problem is to find the “control variable”, u G U, and the “state variable”, (ft G W, 
such that a given cost function, F(u, (ft(u)), defined on IA x W, will be minimized. Here 
(ft satisfies an elliptic PDE which is defined on 0 and will be referred to as the “analysis 
problem” or the “state equation”: 

min ueU F(u,(ft(u)) on T . . 

L(<ft, u) = 0 on Yl 

Note that the control variable is defined on the boundary T, therefore an “optimal boundary 
control” (OBC) problem is considered. 
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2.2 Derivation of the Necessary Conditions for a Minimum 

We apply the adjoint method to the optimal boundary control problem (2.1). The variable 
space is enlarged by adding Lagrange multiplier functions or costate variables denoted by 
A. A Lagrangian is defined to be the sum of the functional and a linear term in the costate 
variables which vanishes as the constraint equation is satisfied; 

A, u ) = F(u, <j>) — (A, L(<f>, u )) . (2.2) 

A perturbation of the Lagrangian with respect to all the variables independently, i.e., state, 
costate and control, results in a variation of the Lagrangian: 

<f> <— 4> + 

A <— A + sA (2.3) 

u u + eu 

with </>, A € L 2 (n), u eU and £ is a small real parameter. The variation of the Lagrange 
function, SE, in the first order approximation in £, is given in the following form: 

SE = — £ (d>, L^(<A, u ) A + F*) - £ (A, L(<j), u)) + £ (u, L* (0, u)A + F u ) (2.4) 

where L * and L* are the adjoint operators of L $ and L u respectively. The requirement 
that the hrst approximation terms vanish results in the necessary condition for a minimum 
which will be referred as the state, the costate, and the design equations: 

State : L(<f), u) = 0 

Costate : u)A + F^{<j), u) = 0 (2.3) 

cojiirol : L’((/>, tt)A + F u (<^>, u) = 0. 

From here on we will use the notation A(u) for the design equation residual, i.e., 

A(u) = -L* u (<l>{u), u)A(u) - F u (<j>(u), u ) (2.6) 

where </>(u) and A(u) in (2.6) are solutions of the state and costate equations. 


2.3 The Sensitivity Gradients 

If the state and costate equations are satisfied, then the variation of the cost function is 
given by (see Eqn.(2.4)): 

8F = -(u,A(u)) r . (2.7) 

This equation implies that the gradient of the functional with respect to the control vari- 
ables is given by — ,4(u): 

V u F(u) = —-4(1*). (2.8) 

Therefore, a perturbation of the control variables with the control residuals multiplied by 
a small parameter, namely u = £«4(u), will result in a reduction of the cost function by 

= -tIMIIi, + 0(£ 2 ). (2.9) 
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2.4 Discretization 


When discretizing the problem it is possible either to derive the necessary conditions for a 
minimum in the continuous formulation and then discretize or to discretize the functional 
together with the state equation and then derive the discrete necessary conditions. In the 
latter case the discrete minimization problem is given by: 

minuh F h (u h , (j> h ) on T h 
L h {<f> h ,u ll ) = 0 onn h . 

As the grid mesh size, h, goes to zero, solutions of both approaches should converge to the 
differential solution. However, for finite mesh size discretization and necessary conditions 
do not necessarily commute. The solutions of both should be within the discretization error 
from the differential solution. For simplicity in this paper we used the first possibility. The 
discrete state, costate and design equations are: 

L h {<j) k ,u h ) = 0 on Sl h 

ii'(^*,«‘)A A + F^,u‘)=0 
£rW*,u‘)A* + F,f(*\ u‘)=0 

We define A h (u h ) similarly to (2.6). 

2.5 A Gradient Descent Algorithm 

The following is a gradient descent minimization algorithm which follows immediately from 
the above. 

1. Start with an initial approximation for the control, Uq . 

2. Solve the state equation for <\> h . 

3. Solve the costate equation for X h . 

4. Compute the amplitude of the perturbation, 0, with a line search, 
and update the control variables: u h <— u h + 0 A h (u' 1 ). 

5. If the residuals of the state, the costate and the control 
equations are greater than some preassigned value, in L 2 norm, 
then goto 2; else stop. 

Note that steps 2, 3 and 5 consist of a global computation over the whole domain. 

The complexity of this algorithm is given by 0(M p N l ), where M is the number of 
control parameters, N is the number of grid points, and p and / are integers which depend 
on the problem and the PDE solver which is used to solve the state and costate equations. 
For example, if a MG solver is used to solve the PDEs then / = 1. 


on Ll h (2.1 1) 

on r\ 
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3 A Multigrid One Shot Minimization Method 

The gradient descent algorithm is applied on a sequence of nested grids, where each coarse 
grid accelerates the convergence rate of its finer grid. On each grid two processes are 
employed: relaxation and coarse grid correction of all the variables, including the control 
variables. On coarse grids the state, the costate and the design equations are restricted 
from the finer grid with the full approximation scheme [3]. 

3.1 Relaxation 

On each level a relaxation is performed on the state, costate and control variables. The 
state and costate equations, which are elliptic PDEs, are relaxed by a Gauss-Seidel or 
damped Jacobi relaxations. The control variables are relaxed by 

u h «- u h + p k f h A h (u k ), (3.1) 

where of and T h are chosen to guarantee good smoothing for the control variables, as 
discussed in Sec. 4, and where A l '(u h ) are the residuals of the design equation. This step 
should be followed by an update of the state and costate solutions. The construction of 
0 h and T h is done so that the boundary data is updated with a high frequency dominated 
quantity. 

In elliptic systems a perturbation of the boundary condition with a Fourier mode e iwx 
has an exponential decaying effect on the interior solution of the form e~ a ^ )y , where y is the 
distance from the boundary and <t(oj) is a positive monotonically increasing function of u> 
for large |w|, [1]. For the Laplace equation the decaying rate is given by e~^ y . Therefore, in 
an MG scheme it is preferable to perturb the boundary condition with only high frequency 
modes relative to the given level. In that case only local relaxations will be needed^in 
order to update the solutions after each optimization step, resulting in an order 0{N~) 
operations for one optimization step. N is the number of interior grid points on a given 
level, and d is the space dimension. On the coarsest grid the relaxation of the control 
variables is given in Sec. 2.5 The PDEs are solved over the whole domain thus taking into 
account the lowest frequencies. In that way the set of grids is complete in the sense that 
all Fourier frequencies are treated at some level. 

3.2 The Coarse Grid Equations 

The restriction of the necessary condition for a minimum to the coarse grid is done with 
the full approximation scheme (see appendix). 

Coarse Grid State Equation 


L H V = /? 

/" = I" si + T t(<f > h ) ° n 

T'K4> h ) = L H iU h - 


(3.2) 



where Iff and Iff are restriction operators which are defined on the interior grid points of 
the domain Q h and which are not necessary identical [3]. 

Coarse Grid Costate Equation 


+ f « = si 

/a ='"/* + »£;(**. A*) 

+ A\u‘)] 

Coarse Grid design equation 


t;"A» + f« = /« 

Si = lift + A‘,«‘) 

A‘, u'“) = A* + F«(/ t V, / 4 "A‘, u‘) 


/» 


L;*A fc + /^(^,A fc ,u fc )] 


on 


-i// 


(3.4) 


where / ff and iff are restriction operators which are defined on the boundary F /l . and 
where the right hand sides /£,/£ and fff are zero on the finest grid. 


3.3 The Coarse Grid Cost Function and Gradient 

It can be shown that the full approximation scheme coarse grid equations, (3. 2-3.4), are the 
necessary conditions for a minimum of the following constrained minimization problem: 

inin^H F H {u H ,<f> H ) - (fff ,^") w - {fff ,u H ) u on V H 

><«) = /« m n", 

where fff, fff and fff are defined in Eqns.(3.2), (3.3) and (3.4). This implies that the 
coarse grid gradient is given by 

Vff„F H = A"(u H ) 

A H (u H ) = fff - (Lf A" + Fff). (3 ' 6) 

Thus the relaxation defined by (3.1), on coarse grids, converges to the solution of the coarse 
grid problem, 

^''(«") = /:-(i; H A» + F«). (3.7) 

3.4 The One Shot Minimization Algorithm 

The problem is solved in one application of an FMG solver. The FMG scheme uses a Vcycle 
scheme in order to solve the problem on each level. The Vcycle is composed of recursive 
applications of a relaxation and coarse grid correction. In the following the relaxation, 
Vcycle and FMG schemes are presented. 
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Relaxation 

A relaxation sweep, 7 Z*\ is defined by the following: 

1. Perform one relaxation of the state equation for (j> h . 

2. Perform one relaxation of the costate equation for A /l . 

3. Update the control variables with the design equation residuals, 
u h «- u h + (3 h T h A h {u h ). 

4. Perform a few local relaxations of the state and costate 
equations in a narrow strip near the boundary. 


Vcycle 

The following is a v 2 ) cycle where V\ and v 2 are integers (in the numerical demon- 

strations we used = 2 and v 2 = 1). The initial grid is the finest, with a mesh size h . 

1. Perform V\ relaxation sweeps, 7 Z h . 

2. Restrict the state, the costate and the design equations 

to the coarse grid (Eqns. (3.2) , (3.3) and (3.4), with H = 2h) . 
Rescale h — ► 2h . 

3. If the coarsest level is not reached goto 1. 

4. Solve the problem with the standard minimization algorithm 
in Sec . 2.5. 

5. Interpolate the coarse grid correction to the finer grid: 

4> h - <t> h + / 2 W - W) ™ 

A /l <- \ h + I% h (\ 2h - ll h \ h ) on n k , 

u h <— u h + / 2 /i( u2/ “ — on ^ k ■ 

Rescale h — * | . 

6. Perform 1/2 relaxation sweeps, H h . 

7. If the finest grid is reached then stop, else goto 5. 

FMG cycle 

The following is a n-FMG(t'i, ^ 2 ) cycle to solve the problem with M grids. The coarsest 
mesh size is denoted by h c . 

1. Start with the coarsest grid, (h = h c ) , and solve the problem with 
the standard minimization algorithm in section 2.5. 

2. Interpolate the solution to a finer grid, rescale h — » | . 

3. Perform n times V h (y\,V 2 ) cycles. 

4. if the finest grid is reached then stop, else goto 2. 
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The computational cost of the cycle is O(N) operations, and it reduces the error of the 
state, costate and control variables by an order of magnitude (in L 2 norm). 

4 Fourier Analysis of the Convergence Rate 

Fourier analysis of the minimization algorithm is described next. The evolution of high 
frequency errors, in the control variables, is considered in half space. Then in a standard 
procedure the problem in half space is reduced to the boundary. A relation between errors 
and residuals of the design equation, on the boundary, is derived. With this relation the 
relaxation and coarse grid correction of the control variable are analyzed. 


4.1 Reduction to a Boundary Problem 

We assume that the state and costate equations are satisfied when the control variables are 
updated. We are interested in the amplification factor of the error in the control variables 
as a result of this process. In the vicinity of the boundary, the non-smooth errors can be 
analyzed using half space geometry. This approximation is valid since in elliptic problems 
non-smooth Fourier modes decay exponentially into the interior (see Sec. 3.1). Consider a 
two-dimensional geometry, where the x-axis is parallel to the boundary and the y- axis is in 
the normal direction. The errors of the state and costate variables satisfy a homogeneous 
equation in the interior at every optimization step, namely 

j> h (9,y = 0 )c ix e l h e- a W y i h dO 

X h {x,y) = JT% A h (9,y = 0 )e ix8 l h e~ s Wyl h dO (4.1) 

u h {x) = u h (0)e ixe / k dd, 

where a(9) and a (6) are determined by the interior state equation: 

^h e ix6/h e -a(e)y/h _ q 

l j h* e ix9/h e -a(e)y/h _ Q (4-2) 

By substituting these expressions into the boundary conditions of the state and costate 
error equations, we obtain relations between <j> h ($,y = 0),A h (9,y = 0) and u h (0 ), which 
are all boundary quantities. Thus, a reduction to a boundary problem has been obtained. 
From the boundary problem we can deduce a relation between the residuals of the design 
equation and the errors in the control variables: 

A h {9) = T h (9)u h (9). (4.3) 

T h (9) is the symbol of the Hessian of the cost function, F, subject to the PDE constraint. 
This symbol determines the smoothing properties of the control variables relaxation as well 
as the effectiveness of the coarse grid correction, as is discussed next. Note that the explicit 
form of the operator, is not known, and in general is a non-local operator. However, 
the computation of its symbol is straightforward. 
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4.2 The Relaxation 


From Eqns.(3.1) and (4.3) it follows that the relation between the errors in the control 
variables before and after the relaxation is given by 

where the relaxation symbol R k (9) is given by 

R h {9) = 1 + {3 ll F h f h {9). (4.5) 

For multigrid proposes it is desirable for R h {9) to have small values in the high-frequency 
range, (| < |0| < 7r). In that case the relaxation will reduce effectively the high-frequency 
errors of the control variables prior to restricting its values to the coarse grid. 

Choice of High Pass Filter 

In some cases the relaxation without the use of a high pass filter (HPF), T h , does 
not smooth the errors effectively for any choice of f3 h . In that case an HPF is introduced 
as a preconditioner of the control residuals. If chosen properly, the symbol T h (9)T h {9) is 
dominated by the high frequencies, and a proper choice of f3 k will result in good smoothing. 
The HPF is particularly effective for problems in which the transformation T h (9) is a 
monotonically decreasing function which has small values in the high frequencies. Without 
the use of a proper HPF, high-frequency oscillatory errors might enter the control variables 
during the computation. 

Evaluation of the optimization step size f} h 

In a multigrid cycle the relaxation should be effective mainly in the high-frequency 
range. The relaxation parameter fi h is chosen so that the maximum of \R h (9)\ in the high 
frequencies will be minimal, that is, 

min max |1 + f} h f h T h (9)\. (4.6) 

flh f <|^|<7T 


(4.4) 


One can show that if the symbol T h {9) does not change sign, then f3 k is given by 


where (j rh T k ) m i n and (T h f h ) 77xax are the minimal and maximal values of J rk (9)f k (9) range 
(§ < |0| < 7 r ) . In most practical problems the symbol P h T h (6) is monotone. Thus (i h is 
given by 

s h = (4 8) 

? h (l)T h (l) + P l {n)T h (ir) 


4.3 Two Level Analysis 

The process for solving the optimal control problem is equivalent to a process of solving the 
equation , T h e h = r h , where e h and r h are the errors and residuals of the design equation, 
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under the assumption that the state and costate equations are satisfied. Using standard 
multigrid arguments the two level convergence matrix is given by (See Appendix) 

M ,i {d) = R h {oy'C h (e)R h {e)^, ( 4 . 9 ) 

where R h {9) is the relaxation matrix 


Rh[e) = ^+^ne) 


1 + (3 h F h T h {0 + tt) 


and C h {6) is the coarse grid correction matrix given by 


W) 

’ r) 


T"( 2»)-' ( /»»(«) , IS (6 + «■))( 


The asymptotic convergence rate is given by the maximum eigenvalue of the matrix M h ($), 
where 0 is in the range 0 < |0| < n. 


5 Numerical Examples 

In this section we demonstrate the performance of the multigrid one shot algorithm and 
apply the analysis developed in Sec. 4, on a series of test problems. The problems are solved 
in a two-dimensional domain which is defined by 

ft = •• o < X < 1 ; 0 < y < l}. 

The constraint is the Poisson equation and the boundary conditions are periodic in 
the z-direction and Dirichlet on the lower boundary, y = 0. The minimization problem is 
defined on the upper boundary, y = 1. 

In subsection 5.1 we solve an optimal control problem of the Dirichlet boundary con- 
dition with four different discretizations. The purpose of this example is to study the 
dependence of the convergence behavior on the choice of discretization. Asymptotic two 
level convergence rates are estimated with Fourier analysis for the different discretization 
schemes. The predicted and the actual convergence rates are compared. 

In subsection 5.2 we solve an optimal control of the Neumann boundary condition for 
two different boundary conditions. The purpose of this example is to show the use of 
the HPF to achieve an efficient smoother for the control variables. The two boundary 
conditions correspond to qualitatively distinct transformations between error and residuals 
of the control equation. The two level analysis is used to determine a proper HPF. 

5.1 The Dirichlet Boundary Control Problem 

Consider the minimization problem is defined by 

min / (^-/*(x)) dx + T) f u 2 dx, (5.1) 

u(x) J y-\ ' on ' Jy=\ 
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where 7 / is a fixed non negative parameter, f*{x) is a given function and where (j> satisfies 
the state equation 

A (f> = / on O 

< <f> = u{x) on y = 1 (5.2) 

k <f> = <t>o on y = 0. 

It is easily verified that the the costate equation is given by 

AA = 0 on S7 

< A + 2(|£ - f*(x)) = 0 on y = 1 (5.3) 

k A = 0 on y = 0 

and the design equation is given by 

4 d\ 

A = -r 2y(f> = 0 on y = 1. (5-4) 

an 


5.1.1 Discretization 


We have used four different discretizations for the minimization problem. For three dis- 
cretizations all unknowns were defined on the vertices the grid lines (referred to as the 
“vertex grid”). The control variables are defined on the intersections of the grid with the 
boundary. In the “cell centered grid” the variables are defined on the centers of the grid 
cells. The grid is extended out of the domain and virtual cell-centered points are defined 
on the neighboring exterior of the domain. A Dirichlet boundary condition is given for the 
average of the variables neighboring the boundary. The control variables are defined on the 
centers of the segments connecting the intersection of the grid with the boundary. Note 
that in the multigrid scheme, the vertices of the grids on different scales are nested while 
in the cell-center case the cells faces are nested. 

In the vertex grid we use three different approximations for the normal derivative on 
the boundary: 

1) A first order approximation for the normal derivative 


VX\ : 


dd 

dm 


— <fc,i 

h 


2) A second order approximation for the normal derivative 


(5.5) 


VX2 : 


8(f) — !<&,i + 2^, 2 — 


(5.6) 


dm h 

3) A use of a virtual point out of the domain, were its value is determined with the 
application of the interior operator on the boundary 


VX3 : 

A cell centered discretization 

CC : 


d(j> _ ^> t ,i - </>«,-! 
dm 2 h 

d^_ _ ~ 

dm h 


(5.7) 


(5.8) 
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5.1.2 Reduction to the Boundary 

In the following we analyze the design equation of the Dirichlet boundary control problem 
in the discrete space. We use a second order finite difference approximation of the Laplacian 
given by 

/ 1 \ 


A* = 


1 


1 

1 -4 1 


V 


1 


(5.9) 


In that case L £ = L and, therefore, cr(9) = cr(9). The term in Eqn.(4.1) satisfies the 
following second order equation (see Eqn.(4.2)) 

+ (-4 + 2 cos 9) + e~ aW = 0. (5.10) 

The Fourier Symbol of the Normal Derivatives 

The normal derivatives, which appear in the design equation, have the following Fourier 
symbols for the different discretizations: 

e-'W - 1 


VXl : 

it* 

II 

VX2 : 

ll 

|co 

VX3 : 

Oil 
3 1 - 

II 

CC: 

d h 

Tn W = 


_I e -M e ) _|_ 2e _tr f e l — | 
h 

p -a(e) _ p a(6) 


(5.11) 

(5.12) 

(5.13) 

(5.14) 
n 

The Fourier Symbol of the Design Equation 

In terms of the normal derivatives the transformation T h (9) (see Eqn.(5.4)) is given by 


2 h 

e-H*> - 


ez 


H«) 


fA w=- 2 i(^wr+4 


dn 


(5.15) 


In this case the calculation of the amplitude of the minimization step, j3 h , given by Eqn.(4.8) 
reduces to 

/3 h = —. 5 -2 • (5.16) 

(£( f )) + (£( 7r )) + 2t i 

In Fig.l the relaxation symbol R h (9) = 1 + 0 ll f h (9) is plotted for the above four 
discretizations. For all four discretizations the relaxation reduces the high frequency errors 
by a factor- smaller than 0.5. 

Fig.2 depicts the maximal eigenvalue, |A| max , of the convergence matrix (4.9) as a 
function of the number of minimization steps, v , on a given level. The factor by which the 
error is reduced as a result of a two level multigrid cycle is bounded by |A| max . It is implied 
by Fig.2 that the cell-centered (CC) and second order vertex (VX2) schemes are expected 
to have a better performance than the other vertex schemes. 
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Figure 1: The symbol of the control variable relaxation for the Dirichlet boundary control 
problem with rj — 0. 
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Figure 2: Two level analysis of asymptotic convergence rates,|A| mox , as a function of the 
number of optimization steps, v, for T} — 0. 
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5.1.3 Convergence Performance 

In the numerical tests the problem (5.1)- (5.2) was solved for the four discretizations (5.5)- 
(5.8). In this problem there was no need to use a high pass filter since the transformation 
T h {9) is dominated by the high frequencies in all four discretizations. The minimization 
step amplitude, given by Eqn.(5.16) was used in the computations. The multigrid one 
shot algorithm was tested using between two and seven levels. The two levels convergence 
is compared with the convergence predicted by the analysis. In all the tests the residuals 
of the state, the costate and the design equations were computed in Z >2 norm. 

In the two levels test the finest grid was composed of 2 7 x 2 7 grid points and the coarsest 
grid was composed of 2 6 x 2 6 grid points. The parameter tj was set to zero. In Fig.3 the two 
level analysis and the actual convergence rates are compared and the similarity between 
them is well apparent. 

In the multilevel test the fine grid was composed of 2 m x 2 m points, with m = 5, 6, 7, and 
the coarsest grid was composed of 2 x 2 grid points. The tests with different choices of m 
were done in order to check if the algorithm is mesh size dependent. All the results in Fig.4 
were done with a cell-centered discretization. Since the case tj = 0 in (5.1) corresponds 
to a trivial problem, the case tj = 1 was tested, although in principle the results should 
not be different. Fig.4 A shows the convergence performance of the analysis problem (5.3). 
Figs.4 B and C show the convergence performance of the optimization problem (5.1) with 
f] = 0 and tj = 1 , respectively. The depicted residuals in 4 B and C are the average of the 
computed state, costate, and design equations residuals. 

In all problems the error was reduced in each Vcycle by an order of magnitude, where 
each Vcycle has a computation complexity of O(N). 
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v v 

Figure 3: Two level convergence rates (C.R.) of the Dirichlet boundary control problem 
as a function of minimization steps on the fine level, v. “TLA” stands for the two level 
analysis prediction and “Numerical” stands for actual convergence rate. The four figures 
correspond to different discretization schemes given by (5.1 1)-(5.14). 
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A. Analysis Problem 



Vcycle Vcycle 

Figure 4: Convergence rates. A. the analysis problem (5.2). B and C the Dirichlet boundary 
control problem, (5.1)-(5.2), with 77 = 0 and 77 = 1 respectively. 


16 





5.2 The Neumann Boundary Control Problem 

In the following examples we compute the Neumann boundary condition on the boundary 
y = \ such that the values of the solution <f> on that boundary will match some given 
function. 


The minimization problem is 

defined by: 




min / ( 6 — 

U(x) Jy= 1 

r{x)) 2 dx, 

(5.17) 

where (f> is satisfying the state equation 




f N> = ! 

on ft, 


i 


on y = 1 

(5.18) 


[ $ = 

on y — 0. 


We study the above problem for 

two different 

right hand sides 

of the boundary condition 


/'i(u) = u 

F 2 {u) = u r . 


(5.19) 


It is easily verified that the the costate equation is 

f AA = 0 on fl 

{ -IS + 2 (* -/’(*)) = 0 on y — \ (5.20) 

[ A = 0 on y — 0. 

The design equation is given on the boundary y = 1 for the two right hand sides, in the 
corresponding order, by 


^4i = — A = 0 

*4.2 ~ A x = 0. 


(5.21) 


Discretization The state and costate variables were discretized on a cell centered grid. 
The control variables are defined on the centers of the segments connecting the intersection 
of the grid with the boundary for the first case: F x {u) — u. In the second case, F 2 (u) = u x , 
the control variables are defined on the intersections of the grid lines with the boundary. 


5.2.1 Analysis 

The symbols of the transformations T h (9) in (4.5) and the proper amplitudes, , calculated 
with Eqn.(4.8) for the different boundary conditions in (5.19) are given by 


fm = 
T 2 h (0) = 


f x 2 yj 6— 2 cos(fl) 

2 1— cos 0 J 


— \/6 — 2 cos 6 ; 


^ h*(s/ 6+V2) . 

^ = TeFTs 


(5.22) 
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One can immediately observe that in the first problem the transformation symbol, T^O), 
is not dominated by the high-frequencies. This means that no choice of can result in 
a relaxation with good smoothing properties, since high-frequency errors will have a low 
weight in the residuals of the design equation. For example, Tj l (ir) = — ^ approaches zero 
with grid refinement; therefore one should expect very slow convergence for high-frequency 
errors. For this case the two level analysis predicts a non-converging scheme, i.e., the 
maximal eigenvalue, |A| maa; , of the convergence matrix (4.9) is greater than one. 

In the second problem there is a high-frequency dominance in the symbol T^iO). There- 
fore one could expect that a few local relaxations near the boundary are needed after each 
perturbation of the control variables. In this case the two level analysis predicts the same 
convergence rate as in the analysis problem, see Fig.6B. 

5.2.2 Use of a High Pass Filter 

In the first test problem, where F\(u) = u, we can perturb the boundary data with, D XX A , 
instead of A , where D xx is a second order tangential derivative. In this case the symbol of 
the design equation will become 

D xx (0)f( l (9) = — 2cos(0), (5.23) 

resulting in a high frequency dominant symbol. The relaxation parameter changes and is 
given by 

P[\D XX ) = -yg + v /g- (5-24) 

A use of a higher order operator such as a fourth order tangential derivative, T = D xxxx , 
results in 

D xxxx (0)Ti(d) = -^6-2008(0X1 - cos(0)) (5.25) 

with 

ft{D xxxx ) = ^ +2v ^- (5-26) 

The two level analysis gives a much better convergence performance for the first choice, see 
Fig.6A. 

5.2.3 Convergence Performance 

The convergence performance of the Neumann boundary control problem is tested for the 
two cases of boundary condition. 

In Fig.7A, the convergence of the optimal control problem (5. 17)-(5. 18), with F(u ) = u, 
is depicted. It is clear from Fig.7A that when using a HPF of the form T = D xx the 
converegnce rate is better than that achieved when using T = D xxxx , as predicted by 
the analysis (Fig.6A). Without using a HPF, {T = /), the algorithm didn’t converge and 
high-frequency oscillatory errors were observed to dominate the solution. 
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A. F(u) = u 


B. F-i{u) = u x 




Figure 5: The symbol of the control variables relaxation, for the Neumann boundary control 
problem (5.17)-(5.18). In Fig. A. two different HPF are used, ((5.23) and (5.25)), with their 
appropriate relaxation parameter (i h , (5.24) and (5.26). 

The two level analysis predicts that if the boundary condition is changed to F(u) = u x 
in Eqn.(5.18) than no HPF 1 is needed for the problem to converge. Fig.7B shows that this 
is indeed the case. 


6 Conclusion 

We have developed a multigrid one shot algorithm to solve the infinite dimensional optimal 
control problem. An analysis, which is based on the reduction of the problem to the 
boundary, was performed both to predict the convergence rate of a two grid algorithm and to 
determine the minimization step. Thus, an expensive line search on every minimization step 
was not required on fine levels. Numerical demonstrations on a series of two dimensional 
test problems were performed. In each test problem the amplitude of the the minimization 
step on fine levels, f) h , and a proper high pass filter (HPF), T h , was determined using 
the analysis. Comparison of the two level convergence rates and two level analysis shows 
agreement within a reasonable error. We find this analysis a simple and powerful tool. In 
each problem, the minimum was reached at a cost of solving the analysis problem just a 
few times, independent of grid size. 
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A. F\(u) = u 
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Figure 6: Two level analysis convergence rates for two choices of boundary conditions in 
the Neumann boundary control problem (see Eqn.(5.19)). T stands for the choice of high 
pass filter. 


A. Optimal Design Problem, F\(u) = u B. Optimal Design Problem, F 2 (u) = u x 



Figure 7: A. Convergence rates for V(2, 1) multigrid cycle of the Neumann boundary control 
problem, with two different HPFs: D xx and D xxxx . Fig. B. No use of HPF is made. Fine 
grid contains 2 7 x ‘2 7 and coarsest 4 grid points. 
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A A Short Review of Multigrid Methods for Solving 
PDEs 

The main idea of multigrid is to solve the problem on a set of nested grids. On each 
grid there are two processes: relaxation and coarse grid correction (CGC). The relaxation 
smooths the errors while the CGC eliminates effectively the smooth errors. In that way 
a very efficient algorithm is achieved which reaches the discretization error in 0(N) op- 
erations, where N is the number of interior grid points. In the following we refer to the 
solution of a discrete elliptic PDE given by 

L h <t> h = f\ (A.l) 

where h is the mesh size. 

Coarse Grid Problem 

For non-linear problems the coarse grid equations are approximated by the full approx- 
imation scheme (FAS): 

L H t H = /,"/'■ + r£(*‘), (A.2) 

where iff and ifa denote the restriction and interpolation operators, respectively, and where 
r£ is defined by 

= (A. 3) 

where If? and if? are not necessary the same [3]. The coarse grid correction (CGC) is given 
by 

$new = fa + Ih{ 4> H ~ Iiftfu)- (A.4) 

Smoothing Properties 

In the infinite space mode analysis the Fourier symbol of the differential operator, 
V\ is denoted by L h (0). For standard MG to work properly one must have h-elliptic 
discretization, i.e. 

ii‘(#)i > c| /<»■ n*i <«• (a.5) 

where L h {6) is the Fourier symbol of L h and where C is a constant. The consequence of 
the condition in (A.5) is that the relaxation will be effective mainly for the high frequency 
errors. For the Jacobi relaxation, the relaxation operator, R h , is related to the difference 
operator, L h , by 

R h = I- pL h , (A. 6) 

where / is the identity operator and ft is a parameter. The smoothing rate of the relaxation 
is determined by 

max |fl‘(«)| < Co < 1, (A.7) 

2 SmS* 

where the constant Co is the predicted smoothing factor (typical value is 0.5). 

Two Level Analysis 
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For the linear case a two level mode analysis can be done to predict the convergence 
rate of the multigrid cycle. The full cycle symbol is give by 

M{9) = k h {6r[I - ll(0)i H (26)-' i“(9)L\e)]R k (dy' , (A. 8) 

where R k stands for the relaxation operator, / is the unit matrix, L H stands for the coarse 
grid operator, and V\ and v-i are integers. 

The convergence rate of the cycle is give by 

/*= sup ||Af(0)||. (A. 9) 

0<|(9|<f 
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